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ABSTRACT 

Axisymmetric stellar wind solutions are presented, obtained by numerically 
solving the ideal magnetohydrodynamic (MHD) equations. Stationary solutions 
are critically analysed using the knowledge of the flux functions. These flux 
functions enter in the general variational principle governing all axisymmetric 
stationary ideal MHD equilibria. 

The magnetized wind solutions for (differentially) rotating stars contain 
both a 'wind' and a 'dead' zone. We illustrate the influence of the magnetic 
field topology on the wind acceleration pattern, by varying the coronal field 
strength and the extent of the dead zone. This is evident from the resulting 
variations in the location and appearance of the critical curves where the wind 
speed equals the slow, Alfven, and fast speed. Larger dead zones cause effective, 
fairly isotropic acceleration to super-Alfvenic velocities as the polar, open field 
lines are forced to fan out rapidly with radial distance. A higher field strength 
moves the Alfven transition outwards. In the ecliptic, the wind outflow is clearly 
modulated by the extent of the dead zone. 

The combined effect of a fast stellar rotation and an equatorial 'dead' 
zone in a bipolar field configuration can lead to efficient thermo-centrifugal 
equatorial winds. Such winds show both a strong poleward collimation and 
some equatorward streamline bending due to significant toroidal field pressure 
at mid-latitudes. We discuss how coronal mass ejections are then simulated on 
top of the transonic outflows. 

Subject headings: methods: numerical — MHD — solar wind — stars: winds, 
outflows 
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1. Introduction 

The solar wind outflow presents a major challenge to numerical modeling since it 
is a fully three-dimensional (3D), time-dependent physical environment, where regions 
of supersonic and subsonic speeds coexist in a tenuous, magnetized plasma. Ulysses 
observations (McComas et al. |1998|) highlighted again that the solar wind about the ecliptic 
plane is fundamentally dynamic in nature, while the fast speed wind across both solar poles 
is on the whole stationary and uniform. Recent SOHO measurements (Hassler et al. |1999D 
demonstrated how the fast wind emanating from coronal holes is rooted to the 'honeycomb' 
structure of the chromospheric magnetic network, making the outflow truly 3D, while the 
daily coronal mass ejections are in essence highly time-varying. Moreover, one really needs 
to study these time- dependent, mult i- dimensional aspects in conjunction with the coronal 
heating puzzle (Holzer & Leer |1997] ). 

Working towards that goal, Wang et al. ( p.998|) recently modeled the solar wind using 
a two-dimensional, time-dependent, magnetohydrodynamic (MHD) description with heat 
and momentum addition as well as thermal conduction. Their magnetic topology shows 
both open (polar) and closed (equatorial) fleld line regions. When heating the closed fleld 
region, a sharp streamer-like cusp forms at its tip as the region continuously expands and 
evaporates. A quasi-stationary wind model results where the emphasis is on reaching a 
qualitative and quantitative agreement with the observed latitudinal variation (reproducing, 
in particular, the sharp transition at roughly ±20° latitude between fast and slow solar 
wind) by tuning the spatial dependence of the artiflcial volumetric heating and momentum 
sources. 

We follow another route towards global solar wind modeling, working our way up 
stepwise from stationary ID to 3D MHD conflgurations. In a pure ideal, stationary, 
axisymmetric MHD approach, numerical simulations can beneflt greatly from analytical 
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theory. This is demonstrated in Ustyugova et al. (|1998|) , where stationary niagneto- 
centrifugally driven winds from rotating accretion disks were calculated numerically and 
critically verified by MHD theory. 



In this paper, we extend the Wang et al. ( |1998|) modeling efforts to 2.5D, by including 
toroidal vector components while remaining axisymmetric. This allows us to explore stellar 
wind regimes where rotation is also important. The magnetic field still has open and 
closed field line regions, but in ideal MHD, the closed field region is a 'dead' zone from 
where no plasma can escape. The unknown coronal heating is avoided by assuming a 
polytropic equation of state and dropping the energy equation all together. The stationary, 
axisymmetric, polytropic MHD models are analysed as in Ustyugova et al. 



In particular, we investigate the effects of (i) having both open and closed field line 
regions in axisymmetric stellar winds; and of (ii) time-dependent perturbations within these 
transonic outfiows. While we still ignore the basic question of why there should be a hot 
corona in the first place, we make significant progress towards fully 3D, dynamic models. 
The advantages of a gradual approach towards such 'final' model were pointed out in 
Keppens & Goedbloed ( |1999a| ). There, we initiated our effort to numerically model stellar 



outflows by gradually relaxing the assumptions inherent in the most well-known solar wind 
model: the isothermal Parker wind (Parker |1958 ). In a sequence of stationary, ID, 1.5D, 



and 2.5D, hydrodynamic and magnetohydrodynamic stellar wind models, all obtained with 
the Versatile Advection Code (VAC, see Toth \m^, |TW7| ; Toth & Keppens |T^; Keppens 



& Toth |1999a| , and [http : / /^aw^ . phys . uu . nl7| toth), we demonstrated that we can now 



routinely calculate axisymmetric magnetized wind solutions for (differentially) rotating 
stars. An important generalization of previous modeling efforts (Sakurai |1985| , |1990|) is that 



the fleld topology can have both open and closed fleld line regions, so we model 'wind' 
and 'dead' zones self-consistently. In essence, our work extends the early model efforts by 
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Pneuman & Kopp ( 1971| ) in (i) going from an isothermal to a polytropic equation of state; 



(ii) allowing for stellar rotation; and (iii) including time-dependent phenomena. While 
we get qualitatively similar solutions for solar-like conditions, we differ entirely in the 
numerical procedure employed and in the way boundary conditions are specified. Keppens 
& Goedbloed ( |1999a|) contained one such MHD model for fairly solar-like parameter values. 
In this paper, we start with a critical examination of this 'reference' model. The obtained 
transonic outflow, accelerating from subslow to superfast speeds, must obey the conservation 
laws predicted by theory, by conserving various physical quantities along streamlines. This 
will be checked in Sect. ^j. Section ^ continues with a physical analysis of the model and 
investigates the influence of the magnetic fleld strength and of the latitudinal extent of 
the 'dead' zone. These parameters have a clear influence on the global wind structure, 
especially evident in the appearance and location of its critical surfaces where the wind 
speed equals the slow, Alfven, and fast magnetosonic speeds. We also present one such 
wind solution for a star which rotates twenty times faster than our sun. Finally, Sect. § 
relaxes the stationarity of the wind pattern, by forcing coronal mass ejections on top of the 
wind pattern. Conclusions are given in Sect. |. 

2. Reference model and conservation laws 

2.1. Solution procedure 

We recall from Keppens & Goedbloed (|1999a|) that we solve the following conservation 
laws for the density p, the momentum vector pv, and the magnetic fleld B: 

fin 

(1) 

(2) 
(3) 



3;+v.(pv)=o. 




^^ + V • [pvv + ptoJ - BB] = 


= Pg, 


(9B 

^- + V ■ (vB - Bv) = 0. 
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Here, ptot = P + \B'^ is the total pressure, I is the identity tensor, g = —{GM^/r^)er is the 
external stellar (mass M*) gravitational field with r indicating radial distance. We assume 
p r^ p"' (dimensionless, we take p = p'^/7), where in this paper we only construct models for 
specified polytropic index 7 = 1.13. This compares to the value 1.05 used in recent work by 
Wu, Guo, & Dryer ( |1997|) and an empirically determined value of 1.46 derived from Helios 1 



data by Totten, Freeman, & Arya ( p.995| ). 

The discretized Eqs. (|l])-® are solved on a radially stretched polar grid in the 
poloidal plane using a Total Variation Diminishing Lax-Friedrich discretization (see e.g. 
Toth & Odstrcil |1996D with Woodward limiting (Collela & Woodward |1984|). Stationary 



{d/dt = 0) solutions are identified when the relative change in the conservative variables 
from subsequent time steps drops below a chosen tolerance (sometimes down to 10^^). 
We explained in Keppens & Goedbloed ( |1999aD how we benefitted greatly from implicit 



time integration (see Toth, Keppens, & Botchev |1998| ; Keppens et al. |1999| ; van der Ploeg, 
Keppens, Toth |1997|) for obtaining axisymmetric {d/dip = 0) hydrodynamic (B = 0) stellar 



outflows characterized by p{R, Z) and \{R,Z), where {R, Z) are Cartesian coordinates 
in the poloidal plane. Denoting the base radius by r*, these hydrodynamic models cover 



r G [l,50]r* and have as escape speed fesc = j2GM^,/r^, = 3.3015cs*, with Cs* the base 
sound speed. They are also characterized by a rotational parameter ( = Vt^^r^^/cs^, = 0.0156 
(if not specifled otherwise), and impose boundary conditions at the base such that (i) 
v^p = (R; and (ii) the poloidal base speed Vp is in accordance with a prescribed radial mass 
flux pvp = /massGr/r^. The valuc of the mass loss rate parameter /mass is taken from a ID 
polytropic, rotating Parker wind valid for the equatorial regions under identical parameter 
values. For ( = 0.0156, we get /mass = 0.01377. We clarify below in which way the values for 
the dimensionless quantities fesc/cs*, C; ^iid /mass relate to the prevailing solar conditions. 

To arrive at a 'reference' MHD wind solution, two more parameters enter the 
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description which quantify the initial field strength and the desired extent of the 'dead' 
zone. A stationary, axisymmetric, magnetized stellar wind is the end result of a time 
stepping process which has the initial density p{R, Z) and toroidal velocity component 
v^{R, Z) from the HD solution with identical 7, fesc, and ( parameters. The poloidal 
velocity is also copied from the HD solution in a polar 'wind' zone where 9 < 6'„ind (upper 
quadrant with ^ = at pole), quantified by its polar angle 6'„ind- The 'dead' zone is 
appropriately initialized by a zero poloidal velocity. The field is initially set to a monopole 
field in the 'wind' zone, where 

Br{R, Z;t = 0)= BoR/r', Bz{R, Z; t = 0) = B^Z/r\ (4) 

where r"^ = R"^ + Z^, coupled to a dipole field in the 'dead' zone with 

Bn{R, Z;t = 0) = Sa^^, BziR, Z;t = 0) = a }^^^ ~^^ . (5) 

The strength of the dipole is taken a^ = i?o/(2cos^wind) to keep the radial field component 
By. continuous at 6^ = ^wind- The initial B^p component is zero throughout. Keppens & 
Goedbloed ( |1999aD took Bq = 3.69 and 6'„ind = 60°, so that the corresponding dead zone 



covered only a ±30° latitudinal band about the stellar equator. For the sun at minimum 
activity, the extent of the coronal hole is typically such that 6'„ind = 30°, so it will be useful 
to vary this parameter in what follows (Sect. |^). We use a resolution of 300 x 40 in the 
full poloidal halfplane, and impose symmetry conditions at both poles, and free outfiow 
at the outer radius 50r* (where all quantities are extrapolated linearly in ghost cells). 
At the stellar base, we similarly extrapolate density and all magnetic field components 
from their initial values, but let these quantities adjust in value while keeping this initial 
gradient in the ghost cells. This implies that the density and the magnetic field at the 
base is determined during the time stepping process to arrive at steady-state. We enforce 
the V ■ B = condition using a projection scheme (Brackbill & Barnes |1980| ), to end up 
with a physically realistic magnetic configuration (despite the 'monopolar' field in the wind 



zone). The stellar boundary condition for the momentum equation allows us to specify a 
differential rotation rate ({9) and latitudinally varying mass flux through /mass (6*) • We set 

PVp = frae.ss{0)er/r^, V^ = C{0)R + B^Vp/ Bp. (6) 

The reference model has a rigid rotation rate according to C = 0.0156, while /mass = 0.01377 
in the wind region and zero in the equatorial dead zone. 



As emphasized in Keppens & Goedbloed ( |1999a| ), our choice of boundary conditions 



is motivated by the variational principle governing all axisymmetric, stationary, ideal 



MHD equilibria (see Sect. 2^). The analytic treatment shows that the algebraic Bernoulli 
equation, together with the cross-field momentum balance, really determine the density 
profile and the magnetic flux function concurrently. In keeping with this formalism, we 
impose a base mass flux and a stellar rotation, and let the density and the magnetic field 
configuration adjust freely at the base. In prescribing the stellar rotation, we exploit 
the freedom available in the variational principle by setting a flux function at the base. 
Noteworthy, the Pneuman & Kopp ( |1971j ) model, as well as many more recent modeling 
efforts for stellar MHD winds, fix the base normal component of the magnetic field together 
with the density. Below, we demonstrate that our calculated meridional density structure 



compares well with recent observations by Gallagher et al. ( |1999|) . 



The values for the dimensionless parameters fesc/cg*, C ^ind Bq (actually the ratio 
of the coronal Alfven speed to c^*) are solar-like in the following sense. At a reference 
radius r^, = 1.25Rq, we take values for the number density No — 10^ cm^"^, temperature 
To = 1.5 X 10^ K, coronal field strength Bo^2G, and rotation rate fi© = 2.998 x 10"*^ s~\ 
For 7 = 1.13 and assuming a mean molecular weight /2 = 0.5, the base sound speed then 
turns out to be c^* = 167.241 km/s, with all dimensionless ratios as used in the reference 
model. Further, the value for the mass loss rate parameter /mass = 0.01377 is then in units 
of 1.06 X 10^^ g/s, so that a split-mo nopole magnetic configuration leads to a realistic mass 
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loss rate M oc 47r/inass — 2.9 x 10~"^^MQyr~-^. Since the reference model has a constant 
mass flux in its wind zone, the presence of the dead zone reduces this value by exactly 
(1 — cos6'„ind) = 1/2. Units enter through the reference radius r*, the base sound speed c^*, 
and the base density p* = NoTripjl (with proton mass nip). 

2.2. Streamfunctions 

The final stationary wind pattern is shown below in Fig. Rl (see also Fig. 5 in Keppens 



& Goedbloed |1999a|) . The physical correctness of this numerically obtained ideal MHD 
solution can be checked as follows. All axisymmetric stationary ideal MHD equilibria are 
derivable from a single variational principle SL = 6 J CdV = with Lagrangian density 



(Goedbloed, Keppens, Lifschitz [1998| ; Keppens & Goedbloed |1999b|) : 



C(M\ ^, V^; R, Z) = -^il - M^) I V^ P -^ + -^, ^. (7) 

To obtain an analytic ideal MHD solution, the minimizing Euler-Lagrange equations need 
to be solved simultaneously for the poloidal flux function ^(-R, Z) and the squared poloidal 
Alfven Mach number M'^{R, Z) = pv^/ B^. Here, Bp = (l/i?)e<^ x V^^. In contrast with the 
translationally symmetric case (Goedbloed & Lifschitz p.997| ; Lifschitz & Goedbloed 1997 ), 



the governing variational principle contains factors i?^, while the profiles Hi and H3 are no 
longer flux functions. In particular, 

n..x'^(^ + ^ + ^), (8) 

H2 = -^^'-'S , (9) 

7-1 

where five flux functions H, fi, S, A, x' enter. These direct integrals of the axisymmetric, 
stationary ideal MHD equations are: 
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• the Bernoulli function (~ energy) 

^(^) = 2"" + ^^1 — —-<+ ^^^^:^' ^^^) 

• the derivative of the stream function x' = dx/dip. Indeed, the poloidal stream 
function x(-R? Z) necessarily obeys xl^); provided that the toroidal component of the 
electric field vanishes. These are immediate checks on the numerical solution, namely 
vrBz = vzBri or the fact that streamlines and field lines in the poloidal plane must 
be parallel (easily seen in Fig. 0). 

• the entropy S, which for our polytropic numerical solutions is constant by construction: 
S ^ 1/7, 

• a quantity related to the angular momentum flux Fam = P'^pRv^p — BpRB^ = pv^A, 

defined as 

A(^) = Rv^ - RBp^^, (12) 

pVp 

• and the derivative of the electric field potential 

fi(^) ^^{v,- ^5,) . (13) 

Various combinations of these flux functions can be made, for instance Goedbloed & 
Lifschitz ( [1997|) used the following flux function (instead of A) 

/(^) = RB^ - Rv/^ = -x'A. (14) 

-Dp 

Figure |l| presents gray-scale contour plots of three streamfunctions / (actually log | / |), 
A, and H, calculated from the reference solution, with its poloidal field lines overlaid. 
Ideally, these contours must match the field line structure exactly: all deviations are due to 
numerical errors. Inspection shows that the agreement is quite satisfactory. A quantitative 
measure of the errors is given in the lower two frames. At bottom left, we plotted the 
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relative deviation of Q from the value enforced at the base ({6) (constant to 0.0156 in this 
model). The solid line marks the 10% deviation, actual values range from [0.0092, 0.0187]. 
It should be noted that this is a very stringent test of the solution, since for the chosen 
parameters, the wind is purely thermally driven, and the stellar rotation is dynamically 
unimportant. The largest deviations are apparent at the rotation axis (symmetry axis) in 
Q, which is not unexpected due to its 1/R dependence. Other inconsistencies are in the 
region which has drastically changed from its initial zero velocity, purely dipole magnetic 
field structure: open field lines coming from the polar regions are now draped around a 
dipolar 'dead' zone of limited radial extent. This dead zone simply corotates with the base 
angular velocity, and has a vanishing poloidal velocity Vp and toroidal field B^. Around 
that zone, the stellar wind traces the open field lines. The final bottom right frame shows 
E^, virtually vanishing everywhere, and only at the very base are values of order O{10~^). 
For completeness, we show the plasma beta f3 = 2p/B'^ = 1 contour which exceeds unity in 
an hourglass pattern that stretches out from the dead zone to large radial distances. 

Overall, the obtained stationary numerical solution passes all criteria for being 
physically acceptable. We expect that most errors disappear when using a higher resolution. 
We already exploited a radial grid accumulated near the stellar surface, necessary for 
resolving the near-surface acceleration. However, we could benefit also from a higher 
resolution in polar angle, now only 40 points for the full half circle, by, for instance, using 
the up-down symmetry. 

3. Extensions of the reference model 

With the accuracy of the numerical solutions confirmed by inspection of the 
streamfunctions, we can start the discussion of the infiuence of the physical parameters Bq, 
^wind, and ( on the global wind structure. First, we present a more detailed analysis of the 
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reference solution itself. 

Fig. 1^ shows the density structure at left, where we plot number density as a function 
of polar angle for three fixed radial distances, namely at the base 1.277?©, at 11. Qi?© and at 
12.7i?0. Keppens & Goedbloed (|1999a| ) already demonstrated the basic effect visible here: 
the equatorial density is higher than the polar density. A recent determination of the (r, 6) 
dependence of the coronal electron densities within IRq < r < 1.2Rq by Gallagher et al. 
( |1999|) concluded that the density falloff is faster in the equatorial region than at the poles, 
and that the equatorial densities within the observed region are a factor of three larger than 
in the polar coronal hole. Their study lists number densities of order lO^cm"'^, as in our 
reference model. Our base density at 1.27Rq (dotted in Fig. ^, left panel) has a distinct 
latitude variation reflecting the combined open-closed field line structure. Interestingly, 
the observations in Gallagher et al. (|1999|) show a similar structure, with quoted values of 



8.3 X lO'^cm"^ at 1.2Rq above the pole, increasing to 1.6 x lO^cm"^ at the same distance 
along the equator. In fact, a dip in the density variation was present due to an active region 
situated above the equator. Qualitatively, we recover this variation at the boundary of the 
dead zone. Again, we stress that the base density is calculated self-consistently, hence not 
imposed as a base boundary condition. The other two radial cuts situated beyond the dead 
zone agree quite well with the conclusions drawn by the observational study. 

The reference wind solution also conforms with some well-known studies in MHD 
wind modeling. Suess & Nerney ( |1973D and Nerney & Suess ( p.975| ) pointed out how 



consistent axisymmetric stellar wind modeling which include magnetic fields and rotation 
automatically lead to a meridional flow away from the equator. At large radial distances, 
the flow profile should be of the form vg oc — sin(2^), with a poleward coUimation of the 
magnetic field. This variation in polar angle is general and independent of the precise base 
field structure. In the middle panel of Fig. ||, we show the latitude dependence of vg at 
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5Oi?0 for the reference model. Note the perfect agreement with the predicted variation. 

Since the solar rotation rate, quantified by the parameter (, is low, the calculated 
wind solution in the poloidal plane should be similar to the one presented by Pneuman & 
Kopp ( |1971| ). They constructed purely poloidal, isothermal and axisymmetric models of 
the solar wind including a helmet streamer (or 'dead' zone). An iterative technique was 
used to solve for the steady coronal expansion, while the density and the radial magnetic 
field were fixed at the base. They enforced a dipolar Br with a strength of IG at the 
poles, half the value we use at the initialization. Their uniform coronal temperature was 
taken to be 1.56 x lO^K, almost identical to our base temperature T^. Their base number 
density was imposed to be 1.847 x lO^cm"^, independent of latitude, and they assumed 
a slightly higher value for the mean molecular weigth, namely /i = 0.608. This leads to 
a base density which is a factor of 2.246 higher than the one used in our model. With 
these differences in mind (together with our polytropic equation of state and the rotational 
effects), we show in the right panel of Fig. ^ the magnetic structure and the location of 
the sonic (where Vp = Cs) and the Alfvenic surface (where Vp = Bp/^) in a manner used 
in the original publication of Pneuman & Kopp (|1971|) , their Fig. 4. In the (1/r, cos6') 
projection, the Alfvenic transition on the equator is at the cusp of the helmet structure 
before the sonic point, while the sonic surface is closer to the solar surface at the poles. The 
qualitative agreement is immediately apparent, although our solution method is completely 
different, most notably in the prescription of the boundary conditions. By calculating the 
base density and magnetic field configuration self-consistently, we generalize the solution 
procedure employed by Pneuman & Kopp ( |1971|) as we gain control of the size of the 
dead zone through our parameter 6'wind- This allows us to study the influence of the base 
topology of the magnetic field on the global wind acceleration pattern in what follows. 

Fig. pi confronts three steady-state wind solutions with our reference model, which 
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differ in the latitudinal extent of the dead zone and/or in the magnetic field strength. With 
^wind = 60° and Bq = 3.69 (corresponding to a 2G base coronal field strength) for the 
reference case A, we increased the latitudinal extent of the dead zone by taking 6'wind = 30° 
in case B, doubled the field strength parameter Bq in model C, and took both Bq = 7.4 
and 6'wind = 30° to arrive at model D. We recall that Bq specifies only the initial field 
strength used in the time-stepping process towards a stationary solution. The final base 
field strength turns out to be of roughly the same magnitude, but differs in its detailed 
latitudinal variation. The changes in the gobal wind pattern are qualified by the resulting 
deformations of the critical surfaces (hourglass curves) where the wind speed equals the 
slow, Alfven, and fast speeds. The plotted region stretches out to ^ 18r^, ~ 22.5Rq. 

By enlarging the dead zone under otherwise identical conditions (from A to B), the 
polar, open field lines are forced to fan out more rapidly with radial distance. As a result, 
the acceleration of the plasma occurs closer to the stellar surface, and the critical curves 
become somewhat more isotropic in polar angle. The Alfven surface moves inwards at the 
poles, and shifts outwards above the now larger dead zone at the equator, approaching a 
circle with an equatorial imprint of the dead zone. If we keep the dead zone small, but 
double the initial field strength Bq (from A to C), the opposite behaviour occurs: the 
critical curves, hence the entire acceleration behaviour of the wind, become much more 
anisotropic. The most pronounced change is an inward shift of the polar slow transition, 
and an outward shift of the Alfven and fast polar transition. This behaviour is in agreement 
with what a Weber-Davis model (Weber & Davis |1967| ) predicts to happen when the field 



strength is increased (note that the Weber-Davis model only applies to the equatorial 
region). When both the field strength and the dead zone are doubled (from A to D), the 
resulting Alfven and fast critical curves are rather isotropic due to the infiuence of the dead 
zone. The polar slow transition is displaced inward while the polar Alfven and fast curve 
are shifted outward, as expected for the higher field. The detailed equatorial behaviour is 
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clearly modulated by the existing dead zone. We note that all wind solutions presented are 
still thermally driven, since the solar-like rotation rate is rather low and the field strengths 
are very modest. The changes are entirely due to reasonable variations in magnetic field 
topology and only a factor of two in field strength. One could tentatively argue that such 
variations occur in the solar wind pattern within its 11-year magnetic cycle. In gray-scale, 
Fig. 1^ shows the absolute value of the toroidal field component | B^p \ (this field changes 
sign across the equator). Note that the stellar rotation has wound up the field lines in a 
zone midway between the poles and the dead zone. For higher rotation rates (see below), 
the associated magnetic pressure build up due to rotation can infiuence the wind pattern 
and cause coUimation (Trussoni, Tsinganos, & Sauty |1997| ). Due to the four- lobe structure, 
one can expect parameter regimes which lead to both poleward coUimation (as in the 
monopole-field models of Sakurai |1985| ) and equatorward streamline bending. 



Figure ^ compares the radial dependence of the poloidal velocity for the four models 
(A, B, C, D) at the pole (left panel) and the equator (middle panel). For comparison, we 
overplotted the same quantities in each panel for a solution with a split- monopole base field 
at the same parameter values. This monopolar field solution is identical in nature to the 
Sakurai ( |1985| , |1990| ) models and was shown in Keppens & Goedbloed ( p.9994 Fig. 4). At 



the pole, a faster acceleration to higher speeds as compared to the reference model results 
from either increasing the field strength or enlarging the dead zone. Moreover, all four 
models show a faster initial acceleration than the corresponding monopolar field model. 
Radio-scattering measurements of the polar solar wind speed (Grail et al. p.996|) indicated 
that the polar wind acceleration is almost complete by IO-R0, much closer than expected. 
Our model calculations show that a fast acceleration can result from modest increases in the 
coronal field strength and dead zone extent (model D has a solar-like dead zone of ±60°). 

The middle panel of Fig. shows the distinct decrease in equatorial wind speed due 
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to the dead zone, when compared to a spht-monopole solution. The equatorial velocities 
are reduced by 10 to 40 km/s, depending on the size of the dead zone and the base field 
strength. Enlarging the dead zone reduces the wind speed significantly (compare A to B, 
and C to D). To a lesser degree, the same effect is true for an increase in coronal field 
strength (compare A to C, and B to D). Keppens & Goedbloed ( |19994 Figure 6) contained 
a polar plot of the velocity and the density at fixed radial distance for the reference model, 
where at least qualitatively, a transition from high density, low speed equatorial wind to 
lower density, high speed polar wind is noticable. As evidenced by Fig. ^ this difference in 
equatorial and polar wind is even more pronounced for larger dead zones. The velocities 
reached are too low for explaining the solar wind speeds - the Weber-Davis wind solution of 
identical parameters reaches 263 km/s at 1 AU. However, this is a well-known shortcoming 



of a polytropic MHD description for modeling the solar wind. Wu et al. (|1999|) therefore 
resort to an ad hoc procedure where the polytropic index is an increasing function of radial 
distance, 7(r), to attain a more realistic 420 km/s wind speed at 1 AU, corresponding to the 
'slow' solar wind. A more quantitative agreement with the observations at these distances 
must await models where we take the energy equation into account and/or model extra 
momentum addition as in Wang et al. ( |1998| ). The equatorial toroidal velocity profile is 
shown at right in Fig. ^. Note that increasing Bq or enlarging the dead zone both negatively 
affect the degree to which the corona corotates with the star. 

Keppens & Goedbloed (|1999a|) also contained a hydrodynamic solution for a much 
faster rotation rate quantified by C = 0.3, or twenty times the solar rotation rate. The 
additional centrifugal acceleration moves the sonic transition closer to the star along the 
equator, and induces an equatorward streamline bending at the base at higher latitudes 
(see also Tsinganos & Sauty |1992|) . One could meaningfully ask what remains of this effect 
when a two-component field structure is present as well. Therefore, we calculate an MHD 
wind for this rotation rate, with Bq = 3.69 and ^wind = 60° as in the reference case. The 
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cor responding mass loss rate parameter (only used in the wind zone) is /mass = 0.01553. 
Fig. ^ displays the wind structure as in Fig. ^, with the gray-scale indicating the logarithm 
of the density pattern. For this solution, we used the up-down symmetry to double the 
resolution in polar angle at the same computational cost. The shape of the critical curves 
has changed dramatically, with a significant outward poleward shift of the Alfven curve, 
and a clear separation between Alfven and fast critical curves - in agreement with a 1.5D 
Weber-Davis prediction. The actual position of the critical surfaces may be influenced by 
an interaction with the outer boundary at 50 r^ in the time-stepping towards a stationary 
solution. In fact, the combined Alfven-fast polar transition has shifted outside the 
computational domain, and the residual could not be decreased to arbitrary small values 
but stagnated at (9(10~^) following this interaction. Within the plotted region of ~ 37Rq, 
the solution is acceptable as explained in Sect. ||. Note how the density structure shows 
an increase towards the equator, causing a very effective thermo-centrifugal acceleration of 
the equatorial wind above the dead zone. The equatorward streamline bending occuring in 
the purely hydrodynamic wind is still important, but now clearly affected by the presence 
of the dead zone. The toroidal magnetic pressure built up by the stellar rotation along the 
mid-latitude open field lines is shaping the wind structure as a whole. In those regions, we 
have pVp/B"^ < 1 together with 2p/B'^ < 1. Thereby, it also leads to streamline bending, 
both poleward as clearly seen in the high latitude field lines, and equatorwards in the 
vicinity of the stellar surface. In this way, the magnetic topology consisting of a dead and a 
wind zone, combined with fast rotation, leads to magnetically dominated collimation along 
the stellar poles, together with magneto-rotational deflections along the equator. The latter 
leads to enhanced densities in the equatorial plane. 

Figure |^ shows the radial dependence at a polar angle 6 = 41.6° of the radial velocity 
Vr, sound speed Cg, radial Alfven speed Ar = Br/^, and azimuthal speed v^. Note 
that the radial velocity reaches up to 400km/s (compare with the ~ 200km/s velocities 
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reached under solar conditions as shown in Fig. ^) , as a result of the additional centrifugal 
acceleration. The rotation plays a significant dynamical role here, in contrast to the 
'solar-like' models discussed earlier and displayed in Fig. |^. The toroidal speed f^ reaches 
above lOOkm/s, a factor of 10 higher than along the ecliptic as shown in Fig. ^ (left panel). 
The corotation obtained by the numerical procedure to find the stationary state can again 
be quantified by the relative error | fi — C I /C- fhe bottom panel of Fig. ^proves that it is 
less than 3 % along that same radial cut. 

4. Triggering coronal mass ejections 

In the process of generating a stationary wind solution, various dynamic phenomena 
take place which may have physical relevance. For instance, the equatorial conic section 
delineated by 6* G [^„ind) tt — ^wind] which was initially static {vp = 0) and dipolar throughout, 
first gets 'invaded' by plasma emanating from the wind zone. The open field structure is 
dragged in towards the equator, and most of the dipolar field is moved out of the domain, 
except for the remaining 'dead' zone. One could qualitatively relate some of these changes 
in the global magnetic topology with observed coronal phenomena. 

In reality however, coronal mass ejections represent major disturbances which happen 
on top of the stationary transonic solar wind. They are associated with sudden, significant 
mass loss and cause violent disruptions of the global field pattern. Most notably, one 
frequently observes the global coronal wind structure to return to its previous stationary 
state, after the passage of the CME. Within the realms of our stellar wind models, we can 
trigger CMEs on top of the outflow pattern, study their motion, and at the same time 
demonstrate that the numerical solutions indeed are stable to such violent perturbations by 
returning to a largely unchanged stationary state. We still restrict ourselves to axisymmetric 
calculations, so the geometry of our 'CME' events is rather artificial. In future work, we 
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intend to model these CMEs in their true 3D setting. 

As background stellar wind, we use a slightly modified model B from the previous 
section. Model B had a large dead zone, 6'wind = 30°, Bq = 3.69, and a rigid rotation with 
( = 0.0156 corresponding to Qq ~ 3 x lO^^s^^. We changed the boundary condition on v^ 
to mimick a 'solar-like' differential rotation, by taking 

CW = CO + C2COS2^ + C4COS4^, (15) 

with (q = 0.0156, (2 = —0.00197, and (4 = —0.00248. This enforces the equator to rotate 
faster than the poles in accord with the observations. As expected for this low rotation 
rate, this has no significant infiuence on the wind acceleration pattern. The coronal mass 
ejection is an equally straigthforward modification of the boundary condition imposed on 
the poloidal momentum equation, namely pVp = fma.ss{0,t)er/r'^ with 

/mass (6*, t) = /wind (6') + 5'CME siu ( ) COS^ ^^^ , (16) 

Vtcme/ V^ ^cme / 

for < t < TcME and ^cme - ocme < d < 9cme + ocme, and otherwise 

/mass(^, ^) = /wmd(^)- (17) 

The wind related mass loss rate /wind(^) contains the polar angle dependence due to the 
dead zone, as before. The extra four parameters control the magnitude of the CME mass 
loss rate gcME, the duration Tcme, and the location ^cme and extent < Ccme < 7r/2 in 
polar angle for the mass ejection. We only present one CME scenario for parameter values 
fi'CME = 2, Tcme = 0.5, ^cme = 60°, and ccme = 30°. Note that the up-down symmetry is 
hereby deliberately broken. This scenario mimicks a mass ejection which detaches from the 
coronal base within 45 minutes and which has an associated mass fiux of about 2 x 10^^ g/s. 
In fact, the total amount of mass lost due to the CME can be evaluated from 

Most = S^fcMETCME— ^ 2 {cOs(^CME — ^CMe) — COs(^CME + OCMe)} • (18) 

^ ~ '^CME 
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For the chosen parameter values, this works out to be M^^^ = ||v^, corresponding to 
0.98 X lO^^g, a typical value for a violent event. 

Figure |^ shows the density difference between the evolving mass ejection and the 
background stellar wind at left, the magnetic field structure (middle), and the toroidal 
velocity component v^p (right) at times t = 1 (Ihr 27' after onset) and t = 3 (4hr 20' after 
onset). The region is plotted up to 15 r^: ~ IS.TSi?©. Although the event is triggered in 
the upper quadrant dead zone only, its violent character also disturbs the overlying open 
field (or wind) zone. The added plasma, trapped in the dead zone, even perturbs the lower 
quadrant wind zone at later times. Note that the CME induces global, abrupt changes in 
the toroidal velocity component. The outermost closed field lines get stretched out radially, 
pulling the dead zone along (see Fig. ^[]^). In ideal MHD calculations, they can never 
detach through reconnection, although numerical diffusion can cause it to happen. We 
observed the outermost field lines of the dead zone to travel outwards without noticable 
reconnection. The overall wind pattern in the first 15 base radii thereby approximately 
returns to its original stationary state, as shown in Figure [^ which gives the solution at 
t = 5 (7hr 13' after onset) and t = 30 (more than 43 hrs after onset). Figure |^ shows 
how a hypothetical spacecraft at 21r* ~ 26.23i?0, close to the ecliptic, would record the 
CME-passage as a sudden increase in density and poloidal velocity which eventually relax 
to their pre-event levels. The event is followed by an increased azimuthal flow regime and 
shows large amplitude oscillations in magnetic field strength and orientation. 

As we assumed axisymmetry, this simulation serves as a crude model for CME-type 
phenomena. Interestingly, axisymmetric numerical simulations of toroidal flux 'belts' 
launched from within the dead zone of a purely meridional, polytropic MHD wind can relate 
favourably to satellite magnetic cloud measurements at lAU (Wu et al. |1999| ). Note that we 
triggered a 'CME' by prescribing a time and space dependent mass flux at the stellar base. 



-21- 



where the density and the magnetic field components could adjust freely. Alternatively, as 
used in studies by Mikic & Linker ( [1994]) , global coronal restructuring can be triggered by 



shearing a coronal arcade. Parameter studies of axisymmetric, but ultimately 3D solutions, 
could investigate the formation and appearance of various MHD shock fronts depending on 
plasma beta, Mach numbers, etc. 

5. Conclusions and outlook 

Continuing our gradual approach towards dynamic stellar wind simulations in three 
dimensions, we studied the infiuence of the magnetic field strength and topology (allowing 
for wind and dead zones), of the stellar rotation, and of sudden mass ejecta on axisymmetric 
MHD winds. 

We demonstrated how reasonable changes in the coronal magnetic field (factor of two 
in field strength and in dead zone extent) infiuence the detailed acceleration behaviour of 
the wind. Larger dead zones cause effective, fairly isotropic acceleration to super-Alfvenic 
velocities since the polar, open field lines are forced to fan out rapidly with radial distance. 
The Alfven transition moves outwards when the coronal field strength increases. The 
equatorial wind outfiow is in these models sensitive to the presence and extent of the dead 
zone, but has, by construction, a vanishing S^ and a /5 > 1 zone from the tip of the dead 
zone to large radial distances. The parameter values for these models are solar-like, hence 
the winds are mostly thermally driven and, in particular, emanated from slowly rotating 
stars. 

For a twenty times faster than solar rotation rate, the wind structure changes 
dramatically, with a clear separation of the Alfven and fast magnetosonic critical curves. At 
these rotation rates, a pure hydrodynamic model predicts equatorward streamline bending 
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from higher latitudes. Our MHD models show how this is now mediated by the dipolar 
dead zone. An equatorial belt of enhanced density stretches from above the dead zone 
outwards, where effective thermo-centrifugally driven outflow occurs. The magnetic field 
structure shows signs of a strong poleward collimation, due to the significant toroidal field 
pressure build up at these spin rates. As pointed out by Tsinganos & Bogovalov ( |1999| ), 
this situation could apply to our own sun in an earlier evolutionary phase. 

It could be of interest to make more quantitative parameter studies of the interplay 
between field topology, rotation rate, etc. in order to apprehend transonic stellar outflows 
driven by combinations of thermal, magnetic, and centrifugal forces. A systematic study of 
the angular momentum loss rates as a function of dead zone extent, magnetic field strength, 
and rotation rate can aid in stellar rotational evolution modeling (Keppens, Charbonneau, 
& McGregor |1995| ; Keppens |1997| ). Specifically, Li (|1999|) pointed out that the present solar 



magnetic breaking rate is consistent with either one of two magnetic topologies: (i) one 
with the standard coronal field strength of ~ 1 G and a small < 2 Rq dead zone; or (ii) 
one with a larger ~ 5 G dipole strength and a sizeable dead zone. When we calculate the 
torque exerted on the star by the magnetized winds A, B, C and D shown in Fig. ^ as 



rwind = 47r/ dOkpr'vR, (19) 

JO 

(noting the axi- and up-down symmetry), we find T^ind — 0.139 x lO'^^dynecm, 
"^wind — 0.062 X lO^Mynecm, r^j^^ ~ 0.246 x lO^Mynecm, and r^^d — 0.123 x lO^Mynecm. 
This confirms Li's result, since a simultaneous doubling of the coronal field strength and 
the dead zone extent (from model A to D) hardly changes the torque magnitude. As could 
be expected, only enlarging the dead zone lowers the breaking efficiency (as pointed out in 



Solanki, Motamen, & Keppens 1997 ), while only raising the field strength leads to faster 
spin-down. Interestingly, a Weber-Davis prediction with governing parameters identical 
to the reference model A (as presented in Keppens & Goedbloed |1999a|) gives a value 
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Twind = 4:7!-pAr\vrA^^*f^A — 2.387 X lO^^dyne cm (all quantities evaluated at the Alfven 
radius r^), one order of magnitude larger! The same conclusion was reached by Priest & 
Pneuman ( [1974|) by estimating the angular momentum loss rate from the purely meridional 
Pneuman & Kopp (|1971|) model. Although the latter model does not include rotation (so 
that A in Eq. (|19D is strictly zero for this model), Priest & Pneuman (|1974|) could estimate 
the torque for a solar rotation rate from the obtained variation of the poloidal Alfven 
radius as a function of latitude (our Fig. Q, right panel). The resulting estimate was only 
15 % of that for a monopole base field. Our exactly evaluated spin-down rates are 2.6 % to 
10.3 % of a split-monopole case. The large difference arises due to the presence of the dead 
zone and the fact that B^ vanishes across the equator for the wind solutions from Fig. ^. 
Indeed, evaluating the torque from Eq. |1^ for the monopolar wind solution from Keppens 
& Goedbloed ( |19994 Figure 4) gives r^md — 2.326 x lO'^^dynecm, in agreement with the 
Weber-Davis estimate. Hence, it should be clear that full MHD modeling is a useful tool to 
further evaluate and constrain different magnetic braking mechanisms. 

We showed how CME events can be simulated on top of these transonic outflows. The 
detailed wind structure is stable to violent mass dumps, even when ejected in the dead zone. 
Note that we restricted ourselves to axisymmetric perturbations, and it will be of interest 
to show whether the axisymmetric solutions are similarly stable to non-axisymmetric 
perturbations (as recently investigated for shocked accretion flows on compact objects in 
Molteni, Toth, & Kuznetsov |1999|) . One could then focus on truly 3D mass ejecta and their 
parametric dependence (possibly allow for direct comparison with LASCO observations of 
coronal mass ejections), or even experiment with unaligned rotation and magnetic axes. 



A 3D time- dependent analytic model by Gibson & Low ( 1998 ) can be used as a further 
check on the numerics. Alternatively, we may decide to zoom in on (3D) details of the 
wind structure at the boundaries of open and closed fleld line regions or about the echptic 
plane, to see whether shear flow driven Kelvin-Helmholtz instabilities (Keppens et al. |1999| ; 
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Keppens & Toth |1999b|) develop in these regions. 



The Versatile Advection Code was developed as part of the project on 'Parallel 
Computational Magneto-Fluid Dynamics', funded by the Dutch Science Foundation (NWO) 
Priority Program on Massively Parallel Computing, and coordinated by JPG. Computer 
time on the Cray C90 was sponsored by the Dutch 'Stichting Nationale Computerfaciliteiten' 
(NCF). We thank Keith MacGregor for suggesting to compare torque magnitudes for 
different models and an anonymous referee for making several useful comments. 
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Fig. 1. — The numerical, ideal MHD stellar wind can be checked to conserve various 
quantities along poloidal streamlines and field lines. For the reference solution, the poloidal 
field lines (solid) must be isolevels in the contourplots of (A) the fiux function / (plotted is 
log I / I); (B) the fiux function A; and (C) the Bernoulli function H. A quantitative error 
estimate is shown in panel (D) where we plot | fi — (^ | /(^, where the solid lines delineate 
dark-shaded regions with deviations > 10 %. Panel (E) shows in a contour plot of E,^, 
that this toroidal electric field component vanishes nearly everywhere. The solid line in (E) 
indicates the /3 = 1 isolevel. 



Fig. 2. — Analysis of the reference solar wind solution. Left panel: Number density from 
pole to pole at three fixed radial distances: 11.9-Rq (solid), 12. Ti?© (dashed), and at the base 
1.27-R0 (dotted and scaled to fit on the figure). Middle panel: the meridional velocity 
component vg as function of polar angle at a fixed SO/?©. Right panel: the magnetic 
field configuration and the poloidal sonic (dashed) and poloidal Alfven (solid) surface as 
in Pneuman & Kopp (1971). 

Fig. 3. — The variation in the detailed wind acceleration pattern due to changes in the 
stellar magnetic field. Poloidal cuts with the star at the centre contain magnetic field lines 
(solid) and poloidal fiow vectors, necessarily parallel to the field lines. The hourglass curves 
indicate the critical curves for slow (dotted), Alfven (solid), and fast (dashed) speeds in the 
wind acceleration. The gray-scale contours indicate the absolute value of the toroidal field 
component | B^p \. Starting from the reference solution in (A), we double the extent of the 
dead zone in (B), raise the field strength to twice its value in (C), and double both the dead 
zone extent and the field strength in (D). 
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Fig. 4. — For the four solutions shown in Fig. ^, we compare at left: the poloidal velocity 
at the pole as a function of radius; middle: the poloidal velocity along the equator; at right: 
the toroidal velocity along the equator. In all three panels, the crosses indicate the same 
quantity for a solution with a split-monopole base field. 

Fig. 5. — A magnetized wind solution for a star rotating 20 times faster than our reference 
'solar' solution. With the star at centre, field lines and poloidal flow vectors are shown at 
right, density contours are given at left, and the critical slow (dotted), Alfven (solid), and 
fast (dashed) curves are shown throughout the poloidal cut, stretching out to 37 Rq. 

Fig. 6. — A purely radial cut through the wind solution from Fig. ^ at a polar angle of 
6 = 41.6°. Top panel shows the radial Alfven speed Ar, sound speed Cg, radial velocity Vr 
and azimuthal velocity v^. The bottom panel confirms the strict corotation achieved: the 
(flux) function Q deviates less than 3 % from its fixed base value (. 

Fig. 7. — A 'coronal mass ejection' simulated on top of an axisymmetric transonic wind. At 
the times indicated (Ihr 27' and 4hr 21' after the onset), we show poloidal cuts of at left: the 
difference in the density pattern between the evolving CME and the original stationary wind 
solution; middle: the poloidal field structure; at right: the toroidal velocity component. 

Fig. 8. — As Fig. |^, times corresponding to 7hr 13' and 43hr past the onset. Note how the 
solution returns to a state almost identical to the original stationary wind solution. 

Fig. 9. — In situ measurement of a CME passage: number density, poloidal velocity, toroidal 
velocity, poloidal field strength, and toroidal field as a function of time at a fixed position of 
26.23i?Q and an angle of 2.25° above the ecliptic. 
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